﻿
inline void klingemberg(double dist, double u0, double rcut, double s1, double s2, double *frep, double *epoten)
{
    double dummy;
    double smed=0.5*(s1+s2);
    if (dist<smed+rcut)
    {


        dummy=u0*exp(-100*(dist-smed));
        ///Para las fuerzas
        *frep=-0.1875*dummy;
        ///Para el potencial
        *epoten=0.00375*dummy;
    }
    else
    {
        *frep=0.;
        *epoten=0.0;
    }
}

inline void magnetic(double dist, double dz, double u0, double *fpartxy, double *fpartz, double *epoten)
{
    double Co, Co2;
    double dist3= dist*dist*dist;
    double dist4= dist3*dist;
    double u30dis=3.0*u0/dist4;
    ///Calculamos la fuerza entre iii y jjj
    Co=dz/dist;
	Co2=Co*Co;

	*fpartxy= u30dis*(5.0*Co2-1.0);
	*fpartz= u30dis*(5.0*Co2-3.0);
	*epoten = u0*(1.0-3*Co2)/(dist3);
}

inline void idealWallInteraction(particula part, double u0, double lboxz,
                                 double cutoff, double *fpart, double *epoten)
{
    ///Distancia de la particula a las paredes superior e inferior
    ///Pared superior: zpu=+lz/2+a
    ///Pared inferior: zpd=-lz/2-a
	double dzu, dzd;
	double cpotHS=0.00375;
	double cpotWall=cpotHS;
	double cforceHS=0.1875;
	double cforceWall=cforceHS;
	double radio = part.radio;
	double dummy;

    dzu=lboxz-(part.pos[2]+0.5*radio);
    dzd=part.pos[2]-0.5*radio;
    //dzu=-part.pos[2]+0.5*lboxz-radio; ///Tiene que ser >=0
    //dzd=-part.pos[2]-0.5*lboxz+radio; ///Tiene que ser <=0

    if(dzu<(0.5*radio+cutoff)) {
        ///Interaccion con pared superior
        dummy=u0*exp(-100.*(dzu-0.5*radio));
        //dummy=u0*exp(-50.*(dzu));
        ///Para las fuerzas
        *fpart=-cforceWall*dummy;
        ///Para el potencial
        *epoten=cpotWall*dummy;
    }
    if(dzd<(radio+cutoff)){

        ///La particula interactua con la pared inferior
        dummy=u0*exp(-50.*(dzd-0.5*radio));
        ///Para las fuerzas
        *fpart=cforceWall*dummy;
        ///Para el potencial
        *epoten=cpotWall*dummy;
    }
}


inline double interaccionSimple(caja *bx)
{
	///Si calcStress = true, se calcula la parte correspondiente al potencial
	///del tensor de esfuerzos microscopicos y se actualiza en el objeto caja
	double massi,massj;
	double fHS, eHS;
	double fMagxy, fMagz, fMagE;
	double dist, dist2, ePoten;
	double fx,fy,fz;
	double virial=0;
	///Tensor de esfuerzos
	double sxy, sxz, syz, szz;
	void stress(caja *bx);
	///Radio de las particulas interactuantes
	double s1=1.0;
	double s2=1.0;
	///Campo magnetico
	double u0 =bx->u0;
	double u0Poly;
	///Ponemos las fuerzas de las particulas de la caja a 0
	//bx->initForces();
	///Cutoff a la interaccion de klingemberg
	double rcutklin=bx->param.rcut;
    massi=massj=1.0;
	///Inicializamos la energia potencial
	ePoten=0.;

	#ifdef CALC_STRESS
	///Inicializamos el tensor de esfuerzos
	sxy=sxz=syz=szz=0.0;
	#endif
    double emag=0.0;

	//Recorremos todos los pares de particulas
	for(int iii=0;iii<bx->npart-1;iii++)
	{
		for(int jjj=iii+1;jjj<bx->npart;jjj++)
		{
			dist2= bx->distancia2(iii,jjj);
			dist=sqrt(dist2);
			#ifdef POLYDISPERSE
			s1=bx->part[iii].radio;
			s2=bx->part[jjj].radio;
			massi=bx->part[iii].mass;
			massj=bx->part[jjj].mass;
			u0Poly=u0*massi*massj;
			#endif

            #ifdef POLYDISPERSE
            klingemberg(dist, u0Poly, rcutklin, s1, s2, &fHS, &eHS);
			magnetic(dist,bx->dr[2], u0Poly, &fMagxy, &fMagz, &fMagE);
			//fMagxy=fMagz=fMagE=0.0;
            #endif

            #ifndef POLYDISPERSE
			klingemberg(dist, u0, rcutklin, s1, s2, &fHS, &eHS);
			magnetic(dist,bx->dr[2], u0, &fMagxy, &fMagz, &fMagE);
			//fMagxy=fMagz=fMagE=0.0;
			#endif
            //fMagxy=fMagz=fMagE=0.0;
			///Calculamos la fuerza entre iii y jjj
            fx= (fMagxy+fHS)*(bx->dr[0])/dist;
            fy= (fMagxy+fHS)*(bx->dr[1])/dist;
            fz= (fMagz+fHS)*(bx->dr[2])/dist;

            #ifdef CALC_STRESS
            sxy+=(bx->dr[0])*fy;
            sxz+=(bx->dr[0])*fz;
            syz+=(bx->dr[1])*fz;
            szz+=(bx->dr[2])*fz;
			#endif

			bx->part[iii].force[0]+=fx;
			bx->part[iii].force[1]+=fy;
			bx->part[iii].force[2]+=fz;

			bx->part[jjj].force[0]-=fx;
			bx->part[jjj].force[1]-=fy;
			bx->part[jjj].force[2]-=fz;

            ///Acualizamos la energia potencial
			ePoten+=eHS+fMagE;
			//eHS=0;
			emag+=fMagE;
			///Actualizamos el virial
			virial-=dist*(fx*bx->dr[0]+fx*bx->dr[1]+fx*bx->dr[2]);

		}
	}
    #ifdef CALC_STRESS
    ///Actualizamos la parte del potencial del tensor de esfuerzos microscopico
    bx->spxy=sxy/bx->vol;
    bx->spxz=sxz/bx->vol;
    bx->spyz=syz/bx->vol;
    bx->spzz=szz/bx->vol;
	#endif
    //cout<<"Em: "<<emag/bx->npart<<"  "<<(ePoten-emag)/bx->npart<<endl;
	bx->epoten=ePoten/bx->npart;
	bx->virial=virial;
	return ePoten;
}

inline double interaccionConfined(caja *bx)
{
    ///Las interacciones consideradas son:
    ///    - Interaccion magnetica entre particulas
    ///    - Interaccion repulsiva particula-particula
    ///    - Interaccion repulsiva particula-pared
	///Si calcStress = true, se calcula la parte correspondiente al potencial
	///del tensor de esfuerzos microscopicos y se actualiza en el objeto caja
    ///En el calculo del tensor de esfuerzos solamente se tienen en cuenta
    ///la interaccion magnetica y la repulsion particula-particula, pero no
    ///la interaccion particula/caja. COMPROBAR!!!
	///Numero de pares interactuantes
	double dist,dist2,ePoten;
	double fx,fy,fz;
	double virial=0;
	///Tensor de esfuerzos
	double sxy, sxz, syz, szz;
	double fcomp, pcomp;
	void stress(caja *bx);
	///Cutoff a la interaccion de klingemberg
	double rcutklin=bx->param.rcut;
	///Radio de las particulas interactuantes
	double s1=1.0;
	double s2=1.0;
	///Campo magnetico
	double u0 = bx->u0;
	///Cutoff interaccion con paredes
	double cutoffwall=bx->param.rcut;
	///Ponemos las fuerzas de las particulas de la caja a 0
	//bx->initForces();
	///Inicializamos la energia potencial
	ePoten=0.;
	///Inicializamos el tensor de esfuerzos
	sxy=sxz=syz=szz=0.0;
	double fHS, eHS;
	double fMagxy, fMagz, fMagE;
	double massi=1.0;
	double massj=1.0;
	double u0Poly;
	///Recorremos todos los pares de particulas
	for(int iii=0;iii<bx->npart-1;iii++){
	    ///Interaccion de la particula iii con las paredes
	    ///Pared superior: zpu=+lz/2+a
	    ///Pared inferior: zpd=-lz/2-a

	    idealWallInteraction(bx->part[iii], u0, bx->lbox[2], cutoffwall, &fcomp, &pcomp);
	    bx->part[iii].force[2]+=fcomp;
		for(int jjj=iii+1;jjj<bx->npart;jjj++){
            #ifdef POLYDISPERSE
			s1=bx->part[iii].radio;
			s2=bx->part[jjj].radio;
			massi=bx->part[iii].mass;
			massj=bx->part[jjj].mass;
			u0Poly=u0*massi*massj;
			#endif

		    ///Distancia con condiciones de contorno periodicas en plano xy pero no en z
			dist2= bx->distancia2(iii,jjj);
			dist=sqrt(dist2);

			#ifdef POLYDISPERSE
            klingemberg(dist, u0Poly, rcutklin, s1, s2, &fHS, &eHS);
			magnetic(dist,bx->dr[2], u0Poly, &fMagxy, &fMagz, &fMagE);
			//fMagxy=fMagz=fMagE=0.0;
            #endif

            #ifndef POLYDISPERSE
			klingemberg(dist, u0, rcutklin, s1, s2, &fHS, &eHS);
			magnetic(dist,bx->dr[2], u0, &fMagxy, &fMagz, &fMagE);
			//fMagxy=fMagz=fMagE=0.0;
			#endif

            ///Calculamos la fuerza entre iii y jjj
            fx= (fMagxy+fHS)*(bx->dr[0])/dist;
            fy= (fMagxy+fHS)*(bx->dr[1])/dist;
            fz= (fMagz+fHS)*(bx->dr[2])/dist;

            #ifdef CALC_STRESS
                sxy+=(bx->dr[0])*fy;
                sxz+=(bx->dr[0])*fz;
                syz+=(bx->dr[1])*fz;
                szz+=(bx->dr[2])*fz;
			#endif

			bx->part[iii].force[0]+=fx;
			bx->part[iii].force[1]+=fy;
			bx->part[iii].force[2]+=fz;

			bx->part[jjj].force[0]-=fx;
			bx->part[jjj].force[1]-=fy;
			bx->part[jjj].force[2]-=fz;

			///Acualizamos la energia potencial
			ePoten+=eHS+fMagE;
			///Actualizamos el virial
			virial-=dist*(fx*bx->dr[0]+fx*bx->dr[1]+fx*bx->dr[2]);

		}
	}
	///Consideramos la interaccion con las paredes de la ultima particula
	///no considerada en el bucle anterior
	 idealWallInteraction(bx->part[bx->npart-1], u0, bx->lbox[2], cutoffwall, &fcomp, &pcomp);

    bx->part[bx->npart-1].force[2]+=fcomp;
	#ifdef CALC_STRESS
		///Actualizamos la parte del potencial del tensor de esfuerzos microscopico
		bx->spxy=sxy/bx->vol;
		bx->spxz=sxz/bx->vol;
		bx->spyz=syz/bx->vol;
		bx->spzz=szz/bx->vol;

	#endif

	bx->epoten=ePoten/bx->npart;
	bx->virial=virial;
	return ePoten;
}

//inline double interaccionConfined(caja *bx)
//{
//    ///Las interacciones consideradas son:
//    ///    - Interaccion magnetica entre particulas
//    ///    - Interaccion repulsiva particula-particula
//    ///    - Interaccion repulsiva particula-pared
//	///Si calcStress = true, se calcula la parte correspondiente al potencial
//	///del tensor de esfuerzos microscopicos y se actualiza en el objeto caja
//    ///En el calculo del tensor de esfuerzos solamente se tienen en cuenta
//    ///la interaccion magnetica y la repulsion particula-particula, pero no
//    ///la interaccion particula/caja. COMPROBAR!!!
//	///Numero de pares interactuantes
//	double dist,dist2,ePoten;
//	double fx,fy,fz;
//	double virial=0;
//	///Tensor de esfuerzos
//	double sxy, sxz, syz, szz;
//	double fcomp, pcomp;
//	void stress(caja *bx);
//	///Cutoff a la interaccion de klingemberg
//	double rcutklin=bx->param.rcut;
//	///Radio de las particulas interactuantes
//	double s1=1.0;
//	double s2=1.0;
//	///Campo magnetico
//	double u0 = bx->u0;
//	///Cutoff interaccion con paredes
//	double cutoffwall=bx->param.rcut;
//	///Ponemos las fuerzas de las particulas de la caja a 0
//	//bx->initForces();
//	///Inicializamos la energia potencial
//	ePoten=0.;
//	///Inicializamos el tensor de esfuerzos
//	sxy=sxz=syz=szz=0.0;
//	double fHS, eHS;
//	double fMagxy, fMagz, fMagE;
//	double massi=1.0;
//	double massj=1.0;
//	double u0Poly;
//	double dr[3];
//	int np=bx->npart;
//	double lx,ly,lz;
//	lx=bx->lbox[0];
//	ly=bx->lbox[1];
//	lz=bx->lbox[2];
//	///Recorremos todos los pares de particulas
//	for(int iii=0;iii<np-1;iii++){
//	    ///Interaccion de la particula iii con las paredes
//	    ///Pared superior: zpu=+lz/2+a
//	    ///Pared inferior: zpd=-lz/2-a
//
//	    idealWallInteraction(bx->part[iii], u0, bx->lbox[2], cutoffwall, &fcomp, &pcomp);
//	    bx->part[iii].force[2]+=fcomp;
//		for(int jjj=iii+1;jjj<np;jjj++){
//            #ifdef POLYDISPERSE
//			s1=bx->part[iii].radio;
//			s2=bx->part[jjj].radio;
//			massi=bx->part[iii].mass;
//			massj=bx->part[jjj].mass;
//			u0Poly=u0*massi*massj;
//			#endif
////
////		    ///Distancia con condiciones de contorno periodicas en plano xy pero no en z
////    dr[0]=bx->part[jjj].pos[0]-bx->part[iii].pos[0];
////    dr[0]-=lx*anint(dr[0]/lx); //Entero + proximo manteniendo signo
////    /*Y*/
////    dr[1]=bx->part[jjj].pos[1]-bx->part[iii].pos[1];
////    dr[1]-=ly*anint(dr[1]/ly); //Entero + proximo manteniendo signo
////    /*Z*/
////    dr[2]=bx->part[jjj].pos[2]-bx->part[iii].pos[2];
////    #ifndef CONFINED
////    dr[2]-=lz*anint(dr[2]/lz); //Entero + proximo manteniendo signo
////    #endif
////
////	dist2=dr[0]*dr[0]+dr[1]*dr[1]+dr[2]*dr[2];
//
//			dist2= bx->distancia2(iii,jjj);
//
//			//dist2=4.;
//			dist=sqrt(dist2);
//
//			#ifdef POLYDISPERSE
//            klingemberg(dist, u0Poly, rcutklin, s1, s2, &fHS, &eHS);
//			magnetic(dist,bx->dr[2], u0Poly, &fMagxy, &fMagz, &fMagE);
//			//fMagxy=fMagz=fMagE=0.0;
//			//fMagxy=fMagz=fMagE=fHS=eHS=0.0;
//            #endif
//
//            #ifndef POLYDISPERSE
//			klingemberg(dist, u0, rcutklin, s1, s2, &fHS, &eHS);
//			magnetic(dist,bx->dr[2], u0, &fMagxy, &fMagz, &fMagE);
//			//fMagxy=fMagz=fMagE=0.0;
//			#endif
//
//            ///Calculamos la fuerza entre iii y jjj
//            fx= (fMagxy+fHS)*(bx->dr[0])/dist;
//            fy= (fMagxy+fHS)*(bx->dr[1])/dist;
//            fz= (fMagz+fHS)*(bx->dr[2])/dist;
//
//            #ifdef CALC_STRESS
//                sxy+=(dr[0])*fy;
//                sxz+=(dr[0])*fz;
//                syz+=(dr[1])*fz;
//                szz+=(dr[2])*fz;
//			#endif
//
//			bx->part[iii].force[0]+=fx;
//			bx->part[iii].force[1]+=fy;
//			bx->part[iii].force[2]+=fz;
//
//			bx->part[jjj].force[0]-=fx;
//			bx->part[jjj].force[1]-=fy;
//			bx->part[jjj].force[2]-=fz;
//
//			///Acualizamos la energia potencial
//			ePoten+=eHS+fMagE;
////			///Actualizamos el virial
////			virial-=dist*(fx*bx->dr[0]+fx*bx->dr[1]+fx*bx->dr[2]);
//
//		}
//	}
//	///Consideramos la interaccion con las paredes de la ultima particula
//	///no considerada en el bucle anterior
//	 idealWallInteraction(bx->part[bx->npart-1], u0, bx->lbox[2], cutoffwall, &fcomp, &pcomp);
//
//    bx->part[np-1].force[2]+=fcomp;
//;
//	#ifdef CALC_STRESS
//		///Actualizamos la parte del potencial del tensor de esfuerzos microscopico
//		bx->spxy=sxy/bx->vol;
//		bx->spxz=sxz/bx->vol;
//		bx->spyz=syz/bx->vol;
//		bx->spzz=szz/bx->vol;
//
//	#endif
//
//	bx->epoten=ePoten/np;
//	bx->virial=virial;
//	return ePoten;
//}

